#include <fstream>

// Local includes
#include "vtk_io_vtk_indip.h"
#include "elem.h"
#include "equation_systems.h"
#include "numeric_vector.h"

using namespace libMesh ;

void VTKIO_NOVTK::cell_connectivity (const Elem* elem, std::vector<unsigned int>& vtk_cell_connectivity)
{
    elem->connectivity(0,VTK,vtk_cell_connectivity);
}

unsigned int VTKIO_NOVTK::cell_offset (const Elem* elem)
{
    std::vector<unsigned int> conn;
    elem->connectivity(0,VTK,conn);
    return conn.size();
}

unsigned int VTKIO_NOVTK::cell_type(const Elem* elem)
{
  unsigned int celltype = 0; // initialize to something to avoid compiler warning
  
  switch(elem->type())
  {
  case EDGE2:				
    celltype = 3; //VTK_LINE
    break;
  case EDGE3:      
    celltype = 21; //VTK_QUADRATIC_EDGE
    break;// 1
  case TRI3:       
    celltype = 5;  //VTK_TRIANGLE
    break;// 3
  case TRI6:       
    celltype = 22;  //VTK_QUADRATIC_TRIANGLE
    break;// 4
  case QUAD4:      
    celltype = 9;  //VTK_QUAD
    break;// 5
  case QUAD8:      
    celltype = 23;  //VTK_QUADRATIC_QUAD
    break;// 6
  case TET4:      
    celltype = 10;  //VTK_TETRA
    break;// 8
  case TET10:      
    celltype = 24; //VTK_QUADRATIC_TETRA
    break;// 9
  case HEX8:    
    celltype = 12;  //VTK_HEXAHEDRON
    break;// 10
  case HEX20:      
    celltype = 25;  //VTK_QUADRATIC_HEXAHEDRON
    break;// 11
  case PRISM6:     
    celltype = 13;   //VTK_WEDGE
    break;// 12
  case PRISM15:   
    celltype = 65;   //VTK_HIGHER_ORDER_WEDGE
    break;// 13
  case PYRAMID5:
    celltype = 14;   //VTK_PYRAMID
    break;// 14
  case QUAD9:      
    celltype = 28;   //VTK_BIQUADRATIC_QUAD
    break;
  case EDGE4:      
  case HEX27:      
  case PRISM18:    
  case INFEDGE2:   
  case INFQUAD4:   
  case INFQUAD6:   
  case INFHEX8:    
  case INFHEX16:   
  case INFHEX18:   
  case INFPRISM6:  
  case INFPRISM12: 
  case NODEELEM:   
  case INVALID_ELEM:
  default:
    {
  	  std::cerr<<"element type "<<elem->type()<<" not implemented"<<std::endl;
  	  libmesh_error();
    }
  }
  return celltype;
}

void VTKIO_NOVTK::write_refined_mesh(bool refined_mesh)
{
   _write_refined_mesh = refined_mesh;
}

void VTKIO_NOVTK::write_equation_systems(const std::string& fname, const EquationSystems& es)
{ 
  const MeshBase& mesh = (MeshBase&)es.get_mesh();
  unsigned int dim = mesh.spatial_dimension();
  MeshBase::const_element_iterator*       p_it  = NULL;
  MeshBase::const_element_iterator* p_end = NULL; 
  MeshBase::const_element_iterator ait  = mesh.active_elements_begin();
  MeshBase::const_element_iterator aend = mesh.active_elements_end(); 
  MeshBase::const_element_iterator lit  = mesh.level_elements_begin(0);
  MeshBase::const_element_iterator lend = mesh.level_elements_end(0); 
  if (_write_refined_mesh) 
  {
     p_it  = & ait;
     p_end = & aend; 
  }
  else 
  {
     p_it  = & lit;
     p_end = & lend; 
  }
  MeshBase::const_element_iterator it = (*p_it);
  MeshBase::const_element_iterator end = (*p_end);

  const unsigned int n_nodes = es.get_mesh().n_nodes();
  std::vector<bool> nodes_bool;
  nodes_bool.assign(n_nodes,false);
  unsigned int n_elems = 0;
  for (it = (*p_it) ; it != end; ++it)
  {
    const Elem *elem  = (*it);
    n_elems++;
    for (unsigned int i = 0; i < elem->n_nodes(); i++)
    {
        nodes_bool[elem->node(i)] = true;
    }
  }
  std::vector<int> nodes_vtk_id;
  nodes_vtk_id.assign(n_nodes, -1);
  unsigned int vtk_id = 0;
  for (unsigned int i = 0; i < nodes_bool.size(); i++)
  {
    if (nodes_bool[i] == true)
    {
       nodes_vtk_id[i] = vtk_id;
       vtk_id++;
    }
  }
  
  // write header
  FILE* pFile=fopen(fname.c_str(),"w");
  fprintf (pFile,"%s","<VTKFile type=\"UnstructuredGrid\"");
  fprintf (pFile,"%s"," version=\"0.1\"");
  fprintf (pFile,"%s\n"," byte_order=\"LittleEndian\">");
  fprintf (pFile,"%s\n","   <UnstructuredGrid>");
  fprintf (pFile,"%s","      <Piece  ");
  fprintf (pFile,"%s%d%s","NumberOfPoints=\"",vtk_id,"\"  ");
  fprintf (pFile,"%s%d%s\n","NumberOfCells=\"",n_elems,"\">");

// write mesh nodes  
  MeshBase::const_node_iterator node_it = mesh.nodes_begin();
  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
  fprintf (pFile,"%s\n","         <Points>");
  fprintf (pFile,"%s%d%s\n","            <DataArray type=\"Float32\" NumberOfComponents=\"",dim,"\" format=\"ascii\">");
  for(node_it=mesh.nodes_begin();node_it!=node_end;++node_it)
  {
    if (nodes_bool[(*node_it)->id()] == true)
    {
      for(unsigned int i = 0; i < dim; i++)
      {
        fprintf (pFile,"%f ",(*node_it)->operator()(i));
      }
    }
  } 
  fprintf (pFile,"%s\n","");
  fprintf (pFile,"%s\n","            </DataArray>");
  fprintf (pFile,"%s\n","         </Points>");
  
////  write solutions 
  std::vector<Number> soln;
  es.build_solution_vector(soln);
  fprintf (pFile,"%s\n","         <PointData>");
  unsigned int var_count = 0;
  const unsigned int n_es_vars = es.n_vars();  
  for(unsigned int i=0;i<es.n_systems();++i) 
  {
	const System& sys = es.get_system(i);
    const unsigned int n_sys_vars = sys.n_vars();
	for(unsigned int j=0;j<n_sys_vars;++j)
    {
      const std::string& varname = sys.variable_name(j);
      fprintf (pFile,"%s%s%s\n","            <DataArray type=\"Float32\" Name=\"",varname.c_str(),"\" NumberOfComponents=\"1\" format=\"ascii\">");
      for(node_it = mesh.nodes_begin();node_it!=node_end;++node_it)
      {
          unsigned int nd_id = (*node_it)->id();
          if (nodes_bool[nd_id] == true)
          {
            fprintf (pFile,"%f ",soln[(nd_id*n_es_vars)+(j+var_count)]);
          }
      }
      fprintf (pFile,"%s\n","");
      fprintf (pFile,"%s\n","            </DataArray>");
    }
    var_count += n_sys_vars;
  }
  fprintf (pFile,"%s\n","         </PointData>");

// write cell data
  fprintf (pFile,"%s\n","         <CellData>");
  fprintf (pFile,"%s\n","         </CellData>");

// write elements connectivity
  fprintf (pFile,"%s\n","         <Cells>");
  fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">");
  for (it = (*p_it) ; it != end; ++it)
  {
    std::vector<unsigned int> vtk_cell_connectivity;  
    const Elem *elem  = (*it);
    this->cell_connectivity(elem,vtk_cell_connectivity);
    for (unsigned int i = 0; i < vtk_cell_connectivity.size(); i++)
    {
       fprintf (pFile, "%d ",nodes_vtk_id[vtk_cell_connectivity[i]]);
    }
  } 
  fprintf (pFile,"%s\n","");
  fprintf (pFile,"%s\n","            </DataArray>");
  fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">");
  unsigned int offset = 0;
  for (it = (*p_it) ; it != end; ++it)
  {
    const Elem *elem  = (*it);
    offset += this->cell_offset(elem);
    fprintf (pFile, "%d ",offset);
  } 
  fprintf (pFile,"%s\n","");
  fprintf (pFile,"%s\n","            </DataArray>");
  fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">");
  for (it = (*p_it) ; it != end; ++it)
  {
    const Elem *elem  = (*it);
    fprintf (pFile, "%d ",this->cell_type(elem));
  } 
  fprintf (pFile,"%s\n","");
  fprintf (pFile,"%s\n","            </DataArray>");
  fprintf (pFile,"%s\n","         </Cells>");
  fprintf (pFile,"%s\n","      </Piece>  ");        
  fprintf (pFile,"%s\n","   </UnstructuredGrid>");
  fprintf (pFile,"%s\n","</VTKFile>");
  fclose (pFile);       
} 
  
void VTKIO_NOVTK::read (const std::string& name)
{
}
  
void VTKIO_NOVTK::write (const std::string& name)
{
}	
  
  
  
